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§ ! ABSTRACT 

^jT ■ We study the time scale T to equipartition in a ID lattice of N masses coupled 

o\ 

by quartic nonlinear (hard) springs (the Fermi-Pasta-Ulam /3 model). We take the 

>V 

initial energy to be either in a single mode 7 or in a package of low frequency 



modes centered at 7 and of width £7, with both 7 and $7 proportional to N. 



o 

These initial conditions both give, for finite energy densities E/N, a scaling in the 



thermodynamic limit (large N) , of a finite time to equipartition which is inversely 



proportional to the central mode frequency times a power of the energy density 
(E/N). A theory of the scaling with (E/N) is presented and compared to the 
numerical results in the range 0.03 < E/N < 0.8. 
PACS numbers: 05.45. +b, 63.20.Ry, 63.10.+a 



I. INTRODUCTION 

In previous work [1,2], the time scale to achieve equipartition in a nonlinear lattice 
of masses coupled by hard springs (FPU-/3 model) was studied with the initial energy 
primarily in a low-frequency mode, of mode number 7. It was shown in [1], numerically 
and theoretically, that energy transfer to high frequency modes is exponentially slow in a 
perturbation (energy) parameter at low energy. The mechanism of a transition to more rapid 
energy transfer is that for energy above a threshold E oc A r_1 , interaction of neighboring low 
frequency modes will lead to a local superperiod beat oscillation, of period Tb oc N 2 /jE 1 , 
that is stochastic. At a sufficiently high energy a resonance with high frequency modes 
leads to a transition to a fast diffusion regime which occurs above a critical energy E = E c , 
independent of N. The equipartition time scale for E > E c was studied numerically in 
[1, 2], with the more detailed numerical investigation of the scaling time to equipartition in 
[2] for 0.2 < E/N < 1.0 giving T oc TbvN. The \JN was interpreted as a size-dependent 
correction to the Tg time scale. In a heuristic treatment of a more general oscillator chain 
we indicated a possible explanation for the yN mode filling factor [3]. If we additionally set 
7 oc iV in Tb, above, then Tb only depends on the total energy per mode (energy density). 
The numerical results in [2] on the other hand, indicated that T, having an additional \/N 
dependence, would become infinite with N at constant energy density, and therefore T 
would not have a finite value in the thermodynamic limit. 

In another study [4], the energy was placed in a low frequency package centered on a 
mode 7 oc N and with an extension £7 oc N. In that work the numerical results indicated 
that the equipartition time was dependent only on the energy density E/N and therefore 
remained finite in the thermodynamic limit, for a finite E/N . However, the measure of 
equipartition used in that study was rather insensitive, so that the exact scalings were 

2 



difficult to obtain. 

Other work has also discussed the question of time to equipartition for various ini- 
tial conditions and nonlinearities [5-10]. In an early investigation, at higher energies, the 
relation between time scales and Lyapunov exponents was investigated [5]. At the higher 
energies investigated in that study, there was primary mode overlap, for which a weak power 
law behavior was observed, with a transition to stochasticity governed by a random matrix 
approximation. The existence of weak chaos at vanishing energy was further explored in 
Ref. [11]. In other studies of more regular systems, such as the discretized integrable sine 
Gordon equation [6], and the FPU-a cubic nonlinearity, which is a two term Taylor series 
approximation to the integrable Toda lattice [8], abrupt transitions were found between 
power law behaviour and "essentially integrable" behaviour. This class of coupled systems 
is well worth further study, but has significantly different properties from the more thor- 
oughly studied FPU-/? system. In another study a (f> A nonlinear chain was introduced to 
allow the transition between power law and exponentially slow behavior to be studied [9]. 
The attention was mainly on the exponentially slow time scales in the low energy (small 
quartic nonlinearity) region, which is not our concern here. In this paper we examine the 
relationship between initial conditions for which the energy is placed primarily in a single 
mode 7, which we call "mechanical" initial conditions, with initial condition for which the 
energy is placed in a finite S'y package of neighboring modes with £7 oc N, which we call 
"physical" initial conditions. We find, numerically, that there is a transient time t oc \/N, 
with mechanical initial conditions, which does not exist or is of minor importance for the 
physical initial conditions. This transient behaviour manifests itself at relatively small N 
(relatively large E/N), explored in [2] and [3] but tends to disappear for large N, thus 
preserving the thermodynamic limit. 
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The FPU-/3 model is a chain of N masses, coupled to nearest neighbors by hard quartic 
nonlinear springs. The Hamiltonian representing the chain is 

H = E f + ^+* ~ *)* + f («H-i - *)*■ C 1 ) 

i=0 

We consider the case of strong springs (3 > 0) and fixed boundaries q = q^+i = 0. The 
constant (3 describing the strength of the anharmonic potential can be scaled to any positive 
value. We vary the energy and fix (3 at the value 0.1 to compare with previous studies of 
the FPU lattice. The equations of motion are integrated using a fourth order symplectic 
integrator. The harmonic part of the Hamiltonian can be put in the form of N independent 
normal modes via the canonical transformation 

N 

qi = ^UijQj j = l,N, (2) 

.7 = 1 

with canonical variables Qj . The columns of the matrix u,,- are the orthonormal eigenvectors 
of the positive definite Hermitian eigenvalue problem for the Q's. The frequencies ujj of the 
normal modes Qj are 



2N + 2 



ujj = 2 sin ( oAr ) . (3) 



The above transformation puts the Hamiltonian (1) into the form 

N 



H = E ( 2^ + ~f Q ^ + (8N + 8) E WiU^wA*, J, k, i)Q t Q 3 Q k Qi, (4) 
where the coefficients G, as calculated in [1], are 

G(i,j, k, £) = Y J B(i + j + k + l), (5) 

p 

where P represents the eight permutations of sign of j, k and I and the function B(x) takes 

the value 1 if the argument is zero, -1 if the argument is ±2(iV + 1), and zero otherwise. 
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II. NUMERICAL CALCULATIONS 

The main numerical tool we use is the calculation of the effective number of modes 
n e ff containing the energy. We use the general formalism of our previous work [1,2] that 
the linear energies Ei = | (Pf + ujfQ'j), i = 1,...,N, are calculated as a function of 
time. The information entropy is given by S = — Y2i=i e i me i> where e, = Ei/ J2i=i Ei are 
the normalized energies. We define the effective number of modes sharing the energy by 
n e ff(t) = exp S [1-3]. We divide n e ft by N, to get a fraction in the range from zero to 
one, which we plot versus time for various values of N and energy densities E/N. We also 
average over 10-20 different realizations of the initial mode phases to minimize fluctuations. 
We take care to distribute the phases of the modes in a random way so that the quartic 
term in (1) does not make the total energy very different from the linear energy. In this way 
one is always close to a set of slightly perturbed linear oscillators as long as (J3E/N)~ 1. 

For "physical" initial conditions we distribute the total energy E uniformly among 
5j = N/16 low-frequency modes ranging from N/64 to 5N/64 with randomly chosen phases, 
for a fixed value of E/N over a range of N, and plot n e ff/N versus time. We compare the 
results to " mechanical" initial conditions in which energy is placed in a fixed number of 
modes where 5~f = 5 with the modes ranging from N/64 to (N/64) + 4. We have worked 
with other values of dj, finding the same qualitative results. 

In Figure 1 we show the evolution of n e ff/N at fixed energy density E/N = 0.1 for 
N ranging from 256 to 4096 for 5j = N/16. The data are seen to lie essentially on a 
universal evolution curve with the correspondence improving for large values of N (error 
bars, where shown, refer to the statistical error computed over the different initial phases; 
otherwise errors are of the size of the symbols). This result is verified numerically for energy 
densities in the range of E/N [0.03,0.8]. In Figure 2 we show the evolution of n e ff/N for 
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the initial E/N = 0.05 and for the same range of N, as in Figure 1, but for £7 = 5 and 
centered around 7 = (N/64) + 2. We see that there is an initial transient with the evolution 
of n e ff/N versus time, different for each N- value, but coalescing later in the evolution at 
longer times and higher values of n e ff/N. If we introduce a factor of y/N to normalize 
the time scale (not shown), we then find that the evolution curves coalesce at early times 
but then diverge. These results are qualitatively consistent with the numerical study in 
[2] (Fig. 5), which indicated an extra volume filling factor, proportional to yN, when the 
energy was placed primarly in a single low frequency mode, typically 7 = 3. In this case 
the primary driving frequency is a beat Qb °^ ^E/N 2 , with 7 = 3, such that the time for 
transferring energy is much longer than in the present situation. This allows the filling of 
the low frequency modes by successive excitations (see ref. [1]) to manifest itself in the 
additional \JN dependence. In Figure 3 we plot n e ff/N as a function of the logarithm of 
time, for a range of E/N values again with "physical" initial conditions (note that n e ff/N 
asymptotically converges to a value which is smaller than 1, due to fluctuations as computed 
in Ref. [6]). The evolution is a monotonically increasing function with an initial transient 
and later an approximately linear increase on the logarithmic time scale. The linear part 
shifts to the left by somewhat less than a decade with every doubling of the energy density, 
which indicates a power law increase of the time scale with (N/E) a having an approximate 
exponent a ~ 3. To be more quantitative we will normalize the time by a function of N/E, 
as in [2] (Fig. 5), to determine if all values of N/E can be fitted on a universal curve, after 
estimating the value of a analytically. 

III. A SCALING ESTIMATE AND NUMERICAL COMPARISON 

In the following we present an approximate theory of Hamiltonian diffusion to explain, 
qualitatively, the power law behavior at low energy densities. We start by assuming that 
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there is an effective number of low-frequency modes 5k that are responsible for stochastically 
transferring energy to the high-frequency modes, leading to equipartition. We assume that 
the initial mode package containing the energy, £7, is such that £7 < 5k, where the size 
of the effective package 5k determines the couplings to high frequency Because, for early 
times, most of the energy is in the low- frequency modes, it is convenient to classify the 
quartic monomial appearing in the sum of (4) depending on how many of the four Q's in it 
belong to low-frequency modes. The largest quartic terms at early times have the four Q's 
of low-frequency modes. These couplings produce deformations of the linear actions of the 
low frequency modes, creating stochastic separatrices for those modes. In [1] it was shown 
that the necessary energy to create separatrices in the low frequency modes has E oc 1/N if 
energy is placed in a single mode. If we place energy in ^7 oc N modes we expect stochastic 
separatrices to be created for E/5j oc 1/N, such that E is independent of N. However, this 
occurs at low energy for 5j/N small. 

Because of the nonlinear couplings among the low- frequency modes, the frequency of 
these modes is corrected by a beat that we approximate in the following way We substitute 
canonical action angle variables Q{ = y/(21i/uji) cos(<pi) and Pi = ^/ ( 2w j ij ) sin (<^>j) into the 
Hamiltonian (4), to obtain 

H = ^2 ^ih + ( 2N 2 ) ^2 G(i,3,k,l)v / uJiUJ j uj k L; l I i I j I k I l ang(ijkl), (6) 

i i,j,k,l 

where ang(ijkl) = cos(^) cos(</>,) cos(^fc) cos(0;). The frequency of mode i is the derivative 
of the Hamiltonian respect to Jj, which evaluates to 



P 






Q i = uj i + [ J 2^ G(i,j, k, tyy/uiVjUJkUiIjIkli ang(ijkl) / \J h. (7) 

• ' j,k,l 

We further assume that there is a rapid spreading over low frequency modes, as observed 
numerically [1], such that we are considering the sum to run over some 5k modes, to be 

7 



determined. After using the selection rule (5), the number of terms in the above sum is 
then of the order of (5k) 2 . We also assume every quartic term in this sum is typically of the 
same size, i.e. with equal energies for all low frequency modes, Wjlj = E/5k. Since these 
terms come with random phases, according to standard gaussian statistics we take the sum 
to be proportional to the square root of the number of terms. With these assumptions, and 
setting u>i = ujj so that (Ij/Ii) 1 / 2 cancels, (7) becomes 

f3E 

ftiRiLJi+LJi— (8) 

with the 5k~ l in the energy per mode cancelling the 5k effective couplings. In (8) and 

below we set N + 1 ~ N (large N) except where it appears in the selection rule. ^,From (3) 

uJi ~ iri/N, for low frequency modes, and taking the beat frequency f2#, between u>i and a 

neighboring mode, to be of the order of the nonlinear frequency shift, we obtain, for any i 

(with i < AT), 

ZIb-- 13 —. (9) 

N N y ' 

The change in the linear energy E^ = \{P 2 + to 2 Q 2 ) of a driving low frequency mode 

i, can be calculated by taking the derivative of (6) with respect to the angle <pi 

~di~ = ~ ( /v+ ) L ° i 5^ G(i,j,hl,h2)^u) i oj j u} hl Lv h2 IiIjIhiIh2sm4> i ang(j,hl,h2), 

(10) 
As in our previous work [1] the notation hi and /i2 explicitly indicates that the energy 
transfer occurs only between low frequency beat oscillation and high frequency mode dif- 
ference oscillations through the Arnold diffusion mechanism. In the above equation, the 
summation is over indices j, hi and /i2 for a given i. The only terms to transfer energy 
to high frequency modes are the ones where j = i, since then the product of the two low 
frequency angles does not have a fast phase associated with it. Additionally, the selection 



rule requires that G(i,j,hl,h2) will be zero unless 

2i + hl + h2 = 2N + 2, (11) 

for which G = — 1. Expanding the dispersion (3) at high frequency and using (11), 

Su h ~ {7r 2 i/2N 2 )(hl - h2). (12) 

In order for the low frequency beat oscillations to transfer energy by Arnold diffusion to 
the high frequency beats we require Q,b ~ Su>h- From (9) and (12) the inequality gives 

-K . j3E > 7T 2 

N % ~W ~ 2iV2 
which reduces to 



i (hi - /i2) Max 



5ft = (hi - h2) Max Z ^E (13) 

To determine the number of terms in the sum in (10) we note that for every i we can take 
h2 arbitrarily from the high frequency package of 5h modes and then hi is calculated from 
(11) with the restriction, from (13) that hi — h2 ~ 2(3E/ir. Writing (11) in the form 

i = N + l-h 2 -£/2 , £ integer, 

such that for i = 1, we have 

/! 2 = AT _ 1, £ = 2 

h 2 = N - 2, I = 4 

up to i = Sh for which 

h 2 = N - Sh/2, £ = Sh. 

Thus we have a decreasing number of couplings with increasing i, with the average number 
of couplings per low frequency mode Sh/4 which scales with f3E as in (13). Substituting 
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this result, together with (13), into (10), then for a single low frequency mode i we obtain 
an estimate for its averaged energy decay 

dE t /2/A J3E „ „ , , 

where from (3) Wj = tti/N. 

Since <5/s low frequency modes, assumed to have energy, couple to Sh high frequency 

modes with 5k = Sh, the cross-couplings imply each high frequency mode is coupled on 

average to Sk/4 low frequency modes. There are phases in the low frequency mode beat 

oscillations and in the high frequency difference oscillations that can affect the Arnold 

diffusion. This has only been studied for exponentially slow diffusion [12]. The effect of 

these phases when more than one driving term exists, for the case of strong Arnold diffusion, 

0,b ~ 8u>h, has not been studied. For lack of evidence we will use the simplest assumption 

that the effect from each low frequency mode is independent. Setting cjj = (3E/N (i = 

Sh/2 = (3E/ir) as an average value in (14), and dividing by E^, we obtain, an average, for 

each mode in the package 

dEi P f(3E\ 2 

with the assumption for scaling purposes that the number of couplings is fixed. Integrating 
(15) in time, with Ei(t) varying from E/Sk at t = to the equipartition value E/N at the 
final time T we get 

ta (»)"@(£) , jf*<<*'- 

Equation (15) only holds, initially, since E decreases in time as the diffusion proceeds. 
However, the change in E is slow compared to the initial build-up of the energy in the 
high frequency modes. Furthermore, we expect that as the energy in the high frequency 
modes increases, other pathways become available for the energy distribution among the 
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modes, to further justify holding the number of couplings constant in the integration. The 

rp 

final step in the approximation is to estimate the value of J Eh(t')dt' at t = T, a time 
of "near-equipartition" . The quantity Eh(t) appears in an integral, so that its exact form 
is not required. For a diffusive process, in which the amplitudes of the modes increase 
with t 1 ' 2 , we might expect the mode energies to increase linearly with t, Eh(t) — ^(E/N), 
such that the time dependence does not depend on TV. This is found to be approximately 
true, numerically, over most of the evolution to near-equipartition. Other forms of the time 
dependence of E^ can also be taken with only small numerical differences. Evaluating the 
integral with the assumption of linear time dependence of Eh (t) we obtain 

T~ 2?F ]n( - V (17) 

We note that the logarithmic factor varies slowly. The numerical coefficient is only a rough 
estimate. Equation (11) exhibits a basic scaling of T a (N/E) 3 . The scaling can be 
checked numerically, by rescaling the time in Fig. 3. This is done in Fig. 4, for five values 
of E/N, giving reasonable confirmation of the rescaling of time with (N/E) 3 . We also 
confirm this scaling by plotting the time to reach n e ft/N = 0.4 against E/N, for all of the 
data, comparing the result to the inverse cubic scaling (dotted line), in Fig. 5. We can also 
compare the magnitude of T in (17) with the numerics. Extrapolating the linear (with log 
time) portion of the E/N = 0.1 curve in Fig. 3 to n e ff/N = 1 we obtain, approximately, 
T ~ 10 . Considering our many approximations, this value is remarkably close to the value 
of T ~ 3 • 10 7 obtained from (17). 

IV. CONCLUSIONS AND DISCUSSION 
We have indicated, numerically, and justified, theoretically, that the FPU-/3 model 
has an appropriate thermodynamic limit. Provided there is sufficient energy in a group 
of low frequency modes that stochastic diffusion to high energy modes occurs on a non- 
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exponentially-slow time scale [1], then the dominant time scale to equipartition is a power 
law (N/E) a . The value of a = 3, estimated from a theoretical scaling argument, was found 
to fit well to the numerical data. 

The numerical results also clarify a result from a previous paper [2] in which a N 1 ' 2 
scaling was numerically found, which would not allow a finite-time thermodynamic limit. 
The resolution of the seeming contradiction is that there is an initial transient which can 
extend over much of the time to equipartition if N is not very large and the initial energy is 
placed in the first few modes (not proportional to N). The existence of a thermodynamic 
limit also agrees with [4], in which the energy was also placed in a mode packet S'j oc N. 
The power of N/E, in that study, numerically fit better to a — 1. The use of a different 
equipartition parameter, less sensitive than n e ff/N, could have led to uncertainty in a, but 
the issue has not been resolved. 

We emphasize that the theory we have developed to explain the scaling, does not predict 
the shape of n e tf{t) which depends on complicated dynamical processes. Furthermore, n e ff 
is related to the evolution of the energy in the individual modes in a very complicated way. 
These dynamics lie beyond a simple mode-averaged theory. We also emphasize that the 
theory depends on having non-exponentially slow stochastic diffusion to high-frequency 
modes, being driven by local mode-mixing stochasticity among low-frequency modes [1]. 
For the approximations to be valid we require that T ^> r where r is the time scale for 
the assumed stochastic process. Since r ~ f^ 1 ~ (N 2 /(3E) 2 , the approximations hold if 
I3E/N < 1. 
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FIGURES 
FIG. 1. rieff/N vs time with E/N = 0.1 and 5-f = N/16, for iV=256, 512, 1024, 2048 and 4096. 

The error bars are the rms variation over (typically) 10 independent trials, which are 

averaged to give final values. 
FIG. 2. n e ff/N vs time with E/N = 0.05 and £7 = 4 for iV=256, 512, 1027, 2048. Error bars 

as in Fig. 1. 
FIG. 3. n e ff/N vs time for TV = 2048 and E/N=0.03, 0.05, 0.1, 0.2, 0.4, and 0.8. 
FIG. 4. rieff/N vs t(/3E/N) 3 . 
FIG. 5. t(n e ff/N = 0.4) vs E/N compared to proportionality t vs (E/N) 3 of dashed line. 
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